Temnos2022Analysis

Required Packages

library(lme4)

## Loading required package: Matrix

library(nlme)

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:lme4':
## 
##     lmList

library(car)

## Loading required package: carData

library(ggplot2)
library(plyr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:nlme':
## 
##     collapse

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(performance)
library(emmeans)
library(devtools)

## Loading required package: usethis

## 
## Attaching package: 'devtools'

## The following object is masked from 'package:emmeans':
## 
##     test

 devtools::install_github("psyteachr/introdataviz")

## Using GitHub PAT from the git credential store.

## Skipping install of 'introdataviz' from a github remote, the SHA1 (0519c980) has not changed since last install.
##   Use `force = TRUE` to force installation

Data

datastore <- read.csv("~/Temnos2022Venom/dissectiondataset.csv")

data<-na.omit(datastore)

summarydata<-read.csv("~/Temnos2022Venom/summaryraiddata2022.csv")

Variables

List:

Data set has both workers and queens. Here are the variables in the dataset:

Date.Dissected: Date ant was dissected MM/DD/YY.

Date.Frozen: Date ant was frozen for dissection. This is approximately 2 weeks after raid trial ended MM/DD/YY.

Date.Frozen: Date colony was collected MM/DD/YY.

Collection.Site: Site at which nest containing ant was collected.

Nest: Number identifier for nest.

Pair.Group: Number identifier for pair group. Each colony collected in the field was divided between 2 nests as a matched pair. Nests with the same Pair group are genetically related (came from the same colony) but one nest in each pair experienced a raiding trial, and one did not.

Raided: Whether or not the nest was raided. Y means yes, N means no.

Worker.Queen : Whether the dissected ant was a worker or a queen

Behavior.Code: If Worker.Queen==” Worker”, this is N for nurse, F for forager, or X for unclassified. If Worker.Queen == “Queen”, this is AQ for alate queen or DQ for dealate queen.

Ant : Ant number (identifier)

Webers.Length.mm : Weber’s length in mm

Venom.Sac.Length.mm : Length of venom sac in mm

Venom.Sac.Width.mm : Width of venom sac in mm

Final.Worker.Count : The number of workers in the nest on the day it was censused/frozen.

Final.Queen.Count.Dealate : The number of queens without wings in the nest on the day it was censused/frozen.

Final.Pupae.Count : The number of pupae in the nest on the day it was censused/frozen.

Final.Larvae.Count : The number of larvae in the nest on the day it was censused/frozen.

Final.Egg.Count : The approximate number of eggs in the nest on the day it was censused/frozen.

Final.Queen.Count.Alate : The number of winged queens in the nest on the day it was censused/frozen.

Final.Male.Count: The number of males in the nest on the day it was censused/frozen.

Add and mutate:

In this chunk, I add a few new variables to the dataset that can be calculated from the initial variables. I also make sure the date is in the correct format. The most important variable created here is Venom.Volume, which is calculated with the formula Venom.Volume \(=\frac{\pi}{6} \times(L \times W^2)\) for each venom sac. I also create a dataset that is only workers (subset of dataset)

data$Raided<-as.factor(data$Raided)
data$Date.Dissected <- as.Date((data$Date.Dissected), format = "%m/%d/%y")
data$Date.Frozen <- as.Date((data$Date.Frozen), format = "%m/%d/%y")
data$Date.Collected <- as.Date((data$Date.Collected), format = "%m/%d/%y")



#In R, first day is automatically Jan 1, 1970. This means that when you put a fixed effect in date format in a model, the results become super difficule to interpret. Therefore we must create a new date variable to include in the model to make results more interpretable. To do this I first convert Frozen date to a numeric, where day 1 is the first day I froze ant nests, August 24, 2022 (since first day in R is 1970 I convert to a numeric and then subtract every day before that, which is 19226 days)
data$Day.Frozen<-as.numeric(data$Date.Frozen)-19226



data$Final.Brood=data$Final.Larvae.Count+data$Final.Pupae.Count
data$Final.Brood.Worker.Ratio=data$Final.Brood/data$Final.Worker.Count

#Scale variables to make results more interpretable based on summary data. Meaning because every ant within a colony has the same final worker count and final pupae count, we need to scale based on the colony level and then give those numbers to the individual level dataset
meanworkers<-mean(summarydata$Final.Worker.Count)
SDworkers<-sd(summarydata$Final.Worker.Count)
data$Worker.Count.Scaled<-(data$Final.Worker.Count-meanworkers)/SDworkers
meanpupae<-mean(summarydata$Final.Pupae.Count)
SDpupae<-sd(summarydata$Final.Pupae.Count)
data$Pupae.Count.Scaled<-(data$Final.Pupae.Count-meanpupae)/SDpupae
meanlarvae<-mean(summarydata$Final.Larvae.Count)
SDlarvae<-sd(summarydata$Final.Larvae.Count)
data$Larvae.Count.Scaled<-(data$Final.Larvae.Count-meanlarvae)/SDlarvae
data$Collection.Site<-as.factor(data$Collection.Site)


#Create a variable for venom volume, which is an ellipsoid using venom sac length and width. Since both length and width are measured in mm, this volume is in mm^3, which is 1:1 with microliters. For the purposes of graphing and modeling this problem, we will multiply by 1000 to convert to nanoliters
data$Venom.Volume<-(pi/6)*data$Venom.Sac.Length.mm*data$Venom.Sac.Width.mm^2*1000






#Create a dataset for just workers
workerdata<-subset(data,Worker.Queen=="Worker")
queendata<-subset(data,Worker.Queen=="Queen")

Does size of ant depend on behavioral caste?

To investigate this we test whether behavior code can predict weber’s length in our data set. We ran two models, one with all ants in the data set and one with just ants from unraided colonies in case ants died during a raid and more ants in the colony switched castes. Regardless, behavioral caste is not a statistically significant predictor of Weber’s length.

 sizemodel<-lme(Webers.Length.mm~Behavior.Code,data=workerdata, random=~1|Nest)
summary(sizemodel)

## Linear mixed-effects model fit by REML
##   Data: workerdata 
##         AIC       BIC   logLik
##   -2062.388 -2040.581 1036.194
## 
## Random effects:
##  Formula: ~1 | Nest
##         (Intercept)   Residual
## StdDev:  0.03227909 0.03746912
## 
## Fixed effects:  Webers.Length.mm ~ Behavior.Code 
##                    Value   Std.Error  DF   t-value p-value
## (Intercept)    0.7968349 0.007648137 552 104.18680  0.0000
## Behavior.CodeN 0.0077771 0.005444028 552   1.42855  0.1537
## Behavior.CodeX 0.0056124 0.005079065 552   1.10500  0.2696
##  Correlation: 
##                (Intr) Bhv.CN
## Behavior.CodeN -0.509       
## Behavior.CodeX -0.537  0.764
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.59849431 -0.59716620  0.03381316  0.64082554  3.01559231 
## 
## Number of Observations: 582
## Number of Groups: 28

plot(sizemodel)

qqnorm(resid(sizemodel))

hist(resid(sizemodel))

 sizemodelunraidedonly<-lme(Webers.Length.mm~Behavior.Code,data=subset(workerdata,Raided=="N"), random=~1|Nest)
summary(sizemodelunraidedonly)

## Linear mixed-effects model fit by REML
##   Data: subset(workerdata, Raided == "N") 
##         AIC       BIC   logLik
##   -1011.311 -992.9269 510.6553
## 
## Random effects:
##  Formula: ~1 | Nest
##         (Intercept)  Residual
## StdDev:  0.03416408 0.0387244
## 
## Fixed effects:  Webers.Length.mm ~ Behavior.Code 
##                    Value   Std.Error  DF  t-value p-value
## (Intercept)    0.7934376 0.011310587 279 70.15000  0.0000
## Behavior.CodeN 0.0055531 0.007816701 279  0.71041  0.4780
## Behavior.CodeX 0.0105755 0.007420036 279  1.42526  0.1552
##  Correlation: 
##                (Intr) Bhv.CN
## Behavior.CodeN -0.503       
## Behavior.CodeX -0.525  0.775
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.32184914 -0.59486995  0.05129671  0.62540029  3.06965088 
## 
## Number of Observations: 295
## Number of Groups: 14

plot(sizemodelunraidedonly)

qqnorm(resid(sizemodelunraidedonly))

hist(resid(sizemodelunraidedonly))

Venom Volume by raided

violinbynest<-ggplot(data, aes(x=as.factor(Nest), y=Venom.Volume))+geom_violin(aes(color=Raided)) + geom_boxplot(width=0.1)

violinbynest

violinraided<-ggplot(subset(data,Worker.Queen=="Worker"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraided

violinraidedforager<-ggplot(subset(workerdata,Behavior.Code=="F"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1) 
violinraidedforager

violinraidednurse<-ggplot(subset(workerdata,Behavior.Code=="N"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraidednurse

violinraidedother<-ggplot(subset(workerdata,Behavior.Code=="X"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraidedother

Analyzing venom volume by raid

#The way Raid.Payoff is currently in dataset is that the raid payoff for the raided colony is also in the corresponding unraided matched pair. This does not make sense to include in the model, so this loop goes through each line and says if the colony didn't experience a raid, the payoff was zero.
for (i in 1:(nrow(workerdata))){
  if (workerdata$Raided[i]=="N"){workerdata$Raid.Payoff.Total[i]=0}
}
workerdata$Final.Pupae.Worker.Ratio<-workerdata$Final.Pupae.Count/workerdata$Final.Worker.Count
workerdata$Final.Brood.And.Eggs<-workerdata$Final.Brood+workerdata$Final.Egg.Count

model<-lme(Venom.Volume~Raided*Behavior.Code+Day.Frozen+ Final.Brood.And.Eggs+Final.Worker.Count+Final.Pupae.Worker.Ratio+I(Webers.Length.mm^3), data=workerdata, random=list(Collection.Site=~1,Pair.Group=~1,Nest=~1))


summary(model)

## Linear mixed-effects model fit by REML
##   Data: workerdata 
##        AIC      BIC    logLik
##   718.4553 783.6661 -344.2276
## 
## Random effects:
##  Formula: ~1 | Collection.Site
##         (Intercept)
## StdDev:  0.02427888
## 
##  Formula: ~1 | Pair.Group %in% Collection.Site
##         (Intercept)
## StdDev:  0.04567579
## 
##  Formula: ~1 | Nest %in% Pair.Group %in% Collection.Site
##         (Intercept)  Residual
## StdDev:  0.07610764 0.4142232
## 
## Fixed effects:  Venom.Volume ~ Raided * Behavior.Code + Day.Frozen + Final.Brood.And.Eggs +      Final.Worker.Count + Final.Pupae.Worker.Ratio + I(Webers.Length.mm^3) 
##                               Value Std.Error  DF   t-value p-value
## (Intercept)               0.1768862 0.1879601 549  0.941084  0.3471
## RaidedY                   0.2914069 0.1051427  10  2.771538  0.0197
## Behavior.CodeN           -0.1870605 0.0828191 549 -2.258664  0.0243
## Behavior.CodeX           -0.1359979 0.0788279 549 -1.725251  0.0850
## Day.Frozen                0.0087471 0.0043938   8  1.990810  0.0817
## Final.Brood.And.Eggs     -0.0008489 0.0011082  10 -0.766080  0.4613
## Final.Worker.Count        0.0017806 0.0025540  10  0.697175  0.5016
## Final.Pupae.Worker.Ratio -1.4126389 1.0496691  10 -1.345794  0.2081
## I(Webers.Length.mm^3)     1.6375319 0.2242739 549  7.301483  0.0000
## RaidedY:Behavior.CodeN   -0.2640293 0.1193151 549 -2.212874  0.0273
## RaidedY:Behavior.CodeX   -0.2406044 0.1111244 549 -2.165180  0.0308
##  Correlation: 
##                          (Intr) RaiddY Bhv.CN Bhv.CX Dy.Frz F.B.A. Fn.W.C
## RaidedY                  -0.282                                          
## Behavior.CodeN           -0.279  0.575                                   
## Behavior.CodeX           -0.262  0.604  0.775                            
## Day.Frozen               -0.532  0.025 -0.014 -0.022                     
## Final.Brood.And.Eggs     -0.328  0.038 -0.033 -0.007  0.731              
## Final.Worker.Count       -0.445  0.014 -0.012 -0.066 -0.045 -0.259       
## Final.Pupae.Worker.Ratio -0.281  0.069  0.051  0.057  0.184  0.134  0.237
## I(Webers.Length.mm^3)    -0.640 -0.032 -0.044 -0.071 -0.003 -0.165  0.228
## RaidedY:Behavior.CodeN    0.199 -0.811 -0.692 -0.538  0.050  0.029  0.015
## RaidedY:Behavior.CodeX    0.201 -0.866 -0.548 -0.706  0.009 -0.008  0.026
##                          F.P.W. I(W.L. RY:B.CN
## RaidedY                                       
## Behavior.CodeN                                
## Behavior.CodeX                                
## Day.Frozen                                    
## Final.Brood.And.Eggs                          
## Final.Worker.Count                            
## Final.Pupae.Worker.Ratio                      
## I(Webers.Length.mm^3)     0.040               
## RaidedY:Behavior.CodeN   -0.049 -0.008        
## RaidedY:Behavior.CodeX   -0.040  0.047  0.767 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.06580499 -0.74442354 -0.06137433  0.53783842  3.76752120 
## 
## Number of Observations: 582
## Number of Groups: 
##                           Collection.Site 
##                                         5 
##           Pair.Group %in% Collection.Site 
##                                        14 
## Nest %in% Pair.Group %in% Collection.Site 
##                                        28

emmeans(model, revpairwise ~ Raided|Behavior.Code)

## $emmeans
## Behavior.Code = F:
##  Raided emmean     SE df lower.CL upper.CL
##  N       1.141 0.0768  4    0.928     1.35
##  Y       1.433 0.0781  4    1.216     1.65
## 
## Behavior.Code = N:
##  Raided emmean     SE df lower.CL upper.CL
##  N       0.954 0.0522  4    0.809     1.10
##  Y       0.982 0.0563  4    0.825     1.14
## 
## Behavior.Code = X:
##  Raided emmean     SE df lower.CL upper.CL
##  N       1.005 0.0448  4    0.881     1.13
##  Y       1.056 0.0447  4    0.932     1.18
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
## Behavior.Code = F:
##  contrast estimate     SE df t.ratio p.value
##  Y - N      0.2914 0.1051 10   2.772  0.0197
## 
## Behavior.Code = N:
##  contrast estimate     SE df t.ratio p.value
##  Y - N      0.0274 0.0703 10   0.389  0.7051
## 
## Behavior.Code = X:
##  contrast estimate     SE df t.ratio p.value
##  Y - N      0.0508 0.0562 10   0.904  0.3871
## 
## Degrees-of-freedom method: containment

joint_tests(model)

##  model term               df1 df2 F.ratio p.value
##  Raided                     1  10   5.545  0.0403
##  Behavior.Code              2 549  14.541  <.0001
##  Day.Frozen                 1   8   3.963  0.0817
##  Final.Brood.And.Eggs       1  10   0.587  0.4613
##  Final.Worker.Count         1  10   0.486  0.5016
##  Final.Pupae.Worker.Ratio   1  10   1.811  0.2081
##  Webers.Length.mm           1 549  53.312  <.0001
##  Raided:Behavior.Code       2 549   2.714  0.0672

plot(model)

qqnorm(resid(model))

hist(resid(model))

#Controlling for body size What is the correct way to control for body size? This chunk runs the same model run in the previous chunk but controls for weber’s length cubed (a measurement of body size) by dividing the explanatory variable by it as opposed to using it as a fixed effect. The results are qualitatively the same and the fit is not as good as in the previous model.

model2<-lme(Venom.Volume/I(Webers.Length.mm^3)~Raided*Behavior.Code+Day.Frozen+ Final.Brood.And.Eggs+Final.Worker.Count+Final.Pupae.Worker.Ratio, data=workerdata, random=list(Collection.Site=~1,Pair.Group=~1,Nest=~1))

summary(model2)

## Linear mixed-effects model fit by REML
##   Data: workerdata 
##        AIC      BIC    logLik
##   1498.201 1559.089 -735.1004
## 
## Random effects:
##  Formula: ~1 | Collection.Site
##         (Intercept)
## StdDev:  0.05681686
## 
##  Formula: ~1 | Pair.Group %in% Collection.Site
##          (Intercept)
## StdDev: 0.0001041914
## 
##  Formula: ~1 | Nest %in% Pair.Group %in% Collection.Site
##         (Intercept)  Residual
## StdDev:   0.1529016 0.8219203
## 
## Fixed effects:  Venom.Volume/I(Webers.Length.mm^3) ~ Raided * Behavior.Code +      Day.Frozen + Final.Brood.And.Eggs + Final.Worker.Count +      Final.Pupae.Worker.Ratio 
##                              Value Std.Error  DF   t-value p-value
## (Intercept)               2.040822 0.2745205 550  7.434132  0.0000
## RaidedY                   0.526210 0.2085407  10  2.523297  0.0302
## Behavior.CodeN           -0.350172 0.1639718 550 -2.135566  0.0332
## Behavior.CodeX           -0.268830 0.1558469 550 -1.724966  0.0851
## Day.Frozen                0.014603 0.0082171   8  1.777152  0.1134
## Final.Brood.And.Eggs     -0.002750 0.0020377  10 -1.349605  0.2069
## Final.Worker.Count        0.005464 0.0044681  10  1.222963  0.2494
## Final.Pupae.Worker.Ratio -3.658864 2.0126765  10 -1.817909  0.0991
## RaidedY:Behavior.CodeN   -0.509998 0.2366036 550 -2.155497  0.0316
## RaidedY:Behavior.CodeX   -0.428976 0.2200864 550 -1.949125  0.0518
##  Correlation: 
##                          (Intr) RaiddY Bhv.CN Bhv.CX Dy.Frz F.B.A. Fn.W.C
## RaidedY                  -0.403                                          
## Behavior.CodeN           -0.416  0.574                                   
## Behavior.CodeX           -0.418  0.604  0.775                            
## Day.Frozen               -0.693  0.021 -0.017 -0.024                     
## Final.Brood.And.Eggs     -0.586  0.029 -0.043 -0.020  0.760              
## Final.Worker.Count       -0.373  0.017  0.000 -0.055 -0.055 -0.226       
## Final.Pupae.Worker.Ratio -0.331  0.066  0.052  0.059  0.181  0.134  0.258
## RaidedY:Behavior.CodeN    0.260 -0.811 -0.693 -0.540  0.056  0.033  0.018
## RaidedY:Behavior.CodeX    0.312 -0.866 -0.547 -0.705  0.012  0.002  0.016
##                          F.P.W. RY:B.CN
## RaidedY                                
## Behavior.CodeN                         
## Behavior.CodeX                         
## Day.Frozen                             
## Final.Brood.And.Eggs                   
## Final.Worker.Count                     
## Final.Pupae.Worker.Ratio               
## RaidedY:Behavior.CodeN   -0.046        
## RaidedY:Behavior.CodeX   -0.040  0.768 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.58930630 -0.70688091 -0.07944495  0.55085331  5.06421913 
## 
## Number of Observations: 582
## Number of Groups: 
##                           Collection.Site 
##                                         5 
##           Pair.Group %in% Collection.Site 
##                                        14 
## Nest %in% Pair.Group %in% Collection.Site 
##                                        28

emmeans(model2, revpairwise ~ Raided|Behavior.Code)

## $emmeans
## Behavior.Code = F:
##  Raided emmean     SE df lower.CL upper.CL
##  N        2.28 0.1510  4     1.86     2.70
##  Y        2.81 0.1537  4     2.38     3.23
## 
## Behavior.Code = N:
##  Raided emmean     SE df lower.CL upper.CL
##  N        1.93 0.1024  4     1.65     2.21
##  Y        1.95 0.1096  4     1.64     2.25
## 
## Behavior.Code = X:
##  Raided emmean     SE df lower.CL upper.CL
##  N        2.01 0.0872  4     1.77     2.25
##  Y        2.11 0.0863  4     1.87     2.35
## 
## Degrees-of-freedom method: containment 
## Results are given on the identity (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
## Behavior.Code = F:
##  contrast estimate    SE df t.ratio p.value
##  Y - N      0.5262 0.209 10   2.523  0.0302
## 
## Behavior.Code = N:
##  contrast estimate    SE df t.ratio p.value
##  Y - N      0.0162 0.139 10   0.116  0.9098
## 
## Behavior.Code = X:
##  contrast estimate    SE df t.ratio p.value
##  Y - N      0.0972 0.112 10   0.871  0.4041
## 
## Note: contrasts are still on the identity scale 
## Degrees-of-freedom method: containment

joint_tests(model2)

##  model term               df1 df2 F.ratio p.value
##  Raided                     1  10   4.203  0.0675
##  Behavior.Code              2 550  13.354  <.0001
##  Day.Frozen                 1   8   3.158  0.1134
##  Final.Brood.And.Eggs       1  10   1.821  0.2069
##  Final.Worker.Count         1  10   1.496  0.2494
##  Final.Pupae.Worker.Ratio   1  10   3.305  0.0991
##  Raided:Behavior.Code       2 550   2.428  0.0892

plot(model2)

qqnorm(resid(model2))

hist(resid(model2))

Visualizing worker data

Behavior Code

Do nurses and foragers have the same amount of venom? Graph it.

venombybehavior<-ggplot(workerdata,aes(Behavior.Code,Venom.Volume,fill=Raided))+
  geom_boxplot()+geom_point(alpha=0.4,position = position_jitterdodge(jitter.width=0.3))+
  xlab("Behavioral Group")+ylab("Venom Volume (nL)")+theme_bw(base_size = 22)+
scale_x_discrete(labels=c('Foragers', 'Nurses', 'Unclassified'))+scale_fill_brewer(palette = "Set2",labels = c("No", "Yes"))

venombybehavior

ggsave('venombybehavior.svg',plot=venombybehavior,dpi=900,width =35, height = 20, units = "cm",device='svg')


venombybehavior<-ggplot(workerdata,aes(Behavior.Code,Venom.Volume,fill=Raided))+
  introdataviz::geom_split_violin(alpha = .4, trim = FALSE) +
  geom_boxplot(width = .2, alpha = .6, fatten = NULL, show.legend = FALSE) +
  stat_summary(fun.data = "mean_se", geom = "pointrange", show.legend = F, 
               position = position_dodge(.175)) +
  xlab("Behavioral Group")+scale_y_continuous(name = "Venom Volume (nL)",
                     breaks = seq(0, 3, 0.5), 
                     limits = c(0, 3.5)) +theme_bw(base_size = 18)+
scale_x_discrete(labels=c('Foragers', 'Nurses', 'Unclassified'))+scale_fill_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes")) +
  theme_minimal()
venombybehavior

## Warning: Removed 233 rows containing missing values or values outside the scale range
## (`geom_split_violin()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).

ggsave('venombybehavior.svg',plot=venombybehavior,dpi=900,width =20, height = 12, units = "cm",device='svg')

## Warning: Removed 233 rows containing missing values or values outside the scale range
## (`geom_split_violin()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).

#Getting warning about NAs, but there are not any
sum(is.na(workerdata$Behavior.Code))

## [1] 0

sum(is.na(workerdata$Venom.Volume))

## [1] 0

sum(is.na(workerdata$Raided))

## [1] 0

#COUNT THE ANTS
raided<-subset(workerdata,Raided=="Y")
unraided<-subset(workerdata,Raided=="N")

nrow(subset(workerdata,Behavior.Code=="N")) #179 nurses

## [1] 179

nrow(subset(raided,Behavior.Code=="N")) #82 raided, 97 unraided

## [1] 82

nrow(subset(unraided,Behavior.Code=="N"))

## [1] 97

nrow(subset(workerdata,Behavior.Code=="F")) #69 foragers

## [1] 69

nrow(subset(raided,Behavior.Code=="F")) #34 raided, 35 unraided

## [1] 34

nrow(subset(unraided,Behavior.Code=="F"))

## [1] 35

nrow(subset(workerdata,Behavior.Code=="X")) #334 unclassified

## [1] 334

nrow(subset(raided,Behavior.Code=="X")) #171 raided, 163 unraided

## [1] 171

nrow(subset(unraided,Behavior.Code=="X"))

## [1] 163

#foragers
ggplot(subset(workerdata, Behavior.Code=="F"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +geom_point()+geom_smooth(method = "lm", fill = NA)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(subset(workerdata, Behavior.Code=="N"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +geom_point()+geom_smooth(method = "lm", fill = NA)

## `geom_smooth()` using formula = 'y ~ x'

Graphical depiction of venom volume by weber’s length

#Graphs by Weber's Length
#overall

weblengthworkers<-ggplot(workerdata, aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +
    geom_point()+geom_smooth(method = "lm", fill = NA)+theme_bw(base_size = 22)+scale_y_continuous(name="Venom volume (nL)", breaks = seq(0, 3, 0.5))+scale_x_continuous(name = ("Weber's length (mm)"),breaks = seq(0.6,1,0.1))+scale_color_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes"))+ggtitle("A")+theme(plot.margin = margin(1, 1, 1, 1, "cm"))
ggsave('weblengthworkers.svg',plot=weblengthworkers,dpi=900,width =25, height = 20, units = "cm",device='svg')

## `geom_smooth()` using formula = 'y ~ x'

weblengthforagers<-ggplot(subset(workerdata,Behavior.Code=="F"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided))+  geom_point()+geom_smooth(method = "lm",se=F)+theme_bw(base_size = 22)+scale_y_continuous(name=expression(paste("Weber's Length (nL)\n",italic("foragers"))), breaks = seq(0, 3, 0.5))+scale_x_continuous(name = "Weber's Length (mm)",breaks = seq(0.6,1,0.1))+scale_color_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes"))+ggtitle("B")+theme(plot.margin = margin(1, 1, 1, 1, "cm"))

ggsave('weblengthforagers.svg',plot=weblengthforagers,dpi=900,width =25, height = 20, units = "cm",device='svg')

## `geom_smooth()` using formula = 'y ~ x'

require(gridExtra)

## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

together<-grid.arrange(weblengthworkers, weblengthforagers, ncol=2)

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

ggsave('weblengthboth.svg',plot=together,dpi=900,width =50, height = 20, units = "cm",device='svg')

This block is for worker venom volume, P1a. First, lets represent graphically.

hist(workerdata$Venom.Volume)

#A bit skewed but mostly pretty normal


venombysize = ggplot(data=workerdata,aes(x=Webers.Length.mm, y=Venom.Volume)) + geom_point() +
  xlab("Weber's Length (mm)")+ylab("Venom Volume")+theme_bw(base_size = 22)+geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE)
#larger ants seem to have more venom, so it's good i'm controlling for that pattern


  
venombydate = ggplot(data=workerdata, aes(Date.Frozen,Venom.Volume)) + 
geom_point() +
geom_smooth(method=lm, se=TRUE ) +
xlab("Date Frozen ")+
ylab("Venom Volume")+
theme_bw(base_size = 22)

Modeling Worker data

by nest

We also need a model to test for overall venom per capita

#Create a column with number of ants dissected in the colony
simplified1<-na.omit(workerdata)
simplified<-simplified1

simplified<-plyr::ddply(simplified, .(Nest), transform, n.Ants.Dissected = length(Ant))
simplified1<-subset(simplified,Behavior.Code=="N")
simplified2<-subset(simplified,Behavior.Code=="F")
#IMPORTANT: THESE DO NOT WORK IF DPLYR HAS BEEN LOADED BEFORE PLYR. BIG ERROR
simplified$Webers.Length.cubed<-(simplified$Webers.Length.mm)^3

simplified<- simplified %>% group_by(Nest, n.Ants.Dissected, Final.Worker.Count,Final.Larvae.Count,Final.Pupae.Count,Final.Queen.Count.Alate,Final.Queen.Count.Dealate,Final.Male.Count, Final.Egg.Count, Day.Frozen,Pair.Group,Raided ) %>%
  summarize(mean.size.cubed=mean(Webers.Length.cubed), mean_venom_volume=mean(Venom.Volume), SizeVar=var(Webers.Length.cubed),ColonyVar=var(Venom.Volume))

## `summarise()` has grouped output by 'Nest', 'n.Ants.Dissected',
## 'Final.Worker.Count', 'Final.Larvae.Count', 'Final.Pupae.Count',
## 'Final.Queen.Count.Alate', 'Final.Queen.Count.Dealate', 'Final.Male.Count',
## 'Final.Egg.Count', 'Day.Frozen', 'Pair.Group'. You can override using the
## `.groups` argument.

nrow(simplified) 

## [1] 28

colonymeanvenommodel<-lme(mean_venom_volume~Raided+mean.size.cubed, data=simplified, random= ~1|Pair.Group)
summary(colonymeanvenommodel)

## Linear mixed-effects model fit by REML
##   Data: simplified 
##         AIC       BIC   logLik
##   -12.44567 -6.351286 11.22283
## 
## Random effects:
##  Formula: ~1 | Pair.Group
##         (Intercept)   Residual
## StdDev:   0.1622871 0.08830954
## 
## Fixed effects:  mean_venom_volume ~ Raided + mean.size.cubed 
##                     Value Std.Error DF  t-value p-value
## (Intercept)     0.3979070 0.3353864 13 1.186414  0.2567
## RaidedY         0.0847105 0.0335927 12 2.521692  0.0268
## mean.size.cubed 1.2019398 0.6372499 12 1.886136  0.0837
##  Correlation: 
##                 (Intr) RaiddY
## RaidedY          0.062       
## mean.size.cubed -0.989 -0.113
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3385189 -0.6473787  0.1547170  0.5135755  1.4247412 
## 
## Number of Observations: 28
## Number of Groups: 14

plot(colonymeanvenommodel)

qqnorm(resid(colonymeanvenommodel))

hist(resid(colonymeanvenommodel))

#A reviewer asked if we need mean size for the colony included in the model. Removing it seems to slightly decrease the p value for Raided but does not change our results qualitatively. It does increase the AIC value, though.
colonymeanvenommodelNOSIZE<-lme(mean_venom_volume~Raided, data=simplified, random= ~1|Pair.Group)
summary(colonymeanvenommodelNOSIZE)

## Linear mixed-effects model fit by REML
##   Data: simplified 
##         AIC       BIC   logLik
##   -10.07063 -5.038247 9.035317
## 
## Random effects:
##  Formula: ~1 | Pair.Group
##         (Intercept)  Residual
## StdDev:   0.1710766 0.0921311
## 
## Fixed effects:  mean_venom_volume ~ Raided 
##                 Value  Std.Error DF   t-value p-value
## (Intercept) 1.0235977 0.05193082 13 19.710794  0.0000
## RaidedY     0.0918653 0.03482228 13  2.638119  0.0205
##  Correlation: 
##         (Intr)
## RaidedY -0.335
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.24495555 -0.40213520 -0.09777063  0.59080947  1.29096471 
## 
## Number of Observations: 28
## Number of Groups: 14

plot(colonymeanvenommodelNOSIZE)

qqnorm(resid(colonymeanvenommodelNOSIZE))

hist(resid(colonymeanvenommodelNOSIZE))

#On average, the venom per capita is higher in raided colonies
jitter <- position_jitter(width = 0.3, height = 0)
venompercapita<-ggplot(simplified, aes(x=Raided, y=mean_venom_volume,fill=Raided))+ylab("Venom per capita (nL)")+geom_boxplot()+geom_point(position=jitter)+theme_bw(base_size = 22)+scale_fill_brewer(palette = "Set2")+theme(legend.position = "None")+scale_x_discrete(labels = c("No", "Yes"))
venompercapita

ggsave('venompercapita2.svg',plot=venompercapita,dpi=900,width =9, height = 20, units = "cm",device='svg')